Analysis of stochastic time series in N dimensions in the presence of strong measurement noise 
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An extension and generalization of a recently presented approach for the analysis of Langevin-type stochastic 
processes in the presence of strong measurement noise is presented. For a stochastic process in A'^ dimensions 
which is superimposed with strong, exponentially correlated, Gaussian distributed, measurement noise it is 
possible to extract the strength and the correlation functions of the noise as well as polynomial approximations 
of the drift and diffusion functions of the underlying process. 
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I. INTRODUCTION 

In the last years there has been significant progress in the 
analysis and characterization of the dynamics of processes un- 
derlying the time series of complex dynamical systems |[l]3. 
If the temporal evolution of a quantity X(i) can be described 
by a Langevin equation, it is possible to extract drift and dif- 
fusion functions of the underlying stochastic process from a 
given timeseries. This can be done because the moments of 
the conditional probability densities of X(<+t) |x(t)=x can be 
related to these functions. 

Since this approach was introduced ||5}i9j it has been suc- 
cessfully carried out in a broad range of fields. For example 
for data from financial markets ifTOl . traffic flow 111], chaotic 
electrical circuits lfT2llT31 . human heart beat ||T4| . climate in- 
dices IIT5I IT6l . turbulent fluid dynamics flT], and for elec- 
troencephalographic data from epilepsy patients ifTsHTOl . 

Real-world data, however, also gives rise to some problems. 
One of them is, that experimental data is only given with a fi- 
nite sampling rate. So methods had to be proposed to deal 
with the effects arising from this fact |20-24|. Another prob- 
lem is the virtually unavoidable measurement noise 13] \24\ - 
l26l . In the presence of measurement noise Y{t) the values of 
X(t) or any of its probability densities are no longer accessi- 
ble, but only X*(t) = X(t) + Y(t) and its density distribu- 
tions. 

Recently an approach has been presented which allows the 
estimation of drift and diffusion functions in the presence 
of strong, delta-correlated, Gaussian noise ETl l28l . Start- 
ing with initial estimates for the noise strength and the drift 
and diffusion functions a functional of these unknowns is it- 
eratively minimized. An alternative approach, which is able 
to deal also with strong, exponentially correlated, Gaussian 
noise, has been presented in ||29l . 

The aim of this paper is the formulation of this later ap- 
proach in N dimensions and also its generalization. The basic 
idea stays the same. Instead of looking at the conditional mo- 
ments in the first place, the joint probability density p(x, x', t) 
of pairs (X(t), X(t+r)) is looked at. If the measurement noise 
is independent of X(i), then (X(i), X(i+r)) and (Y(t), Y(t+ 
t)) are independent random variables. Hence the joint prob- 
ability density p*(x,x', t) of their sum (X*(t),X*(t + T) is 
given by the convolution of p and py, where py (x, x', r) is 
the joint probabiHty density of (Y(t), Y(t+T)). 



The noise is assumed to be Gaussian and the Gauss func- 
tion has some special algebraic properties. This allows to ex- 
press the moments of p* in terms of the moments of p and 
of the noise parameters. The obtained relations can then be 
used to extract the noise parameters. Furthermore, by the use 
of integral transformes (the Fourier transform used in |29| is 
a special case hereof), it is possible to extract polynomial ap- 
proximations of the drift and diffusion functions using purely 
algebraic relations between quantities that can be calculated 
directly from a given, noisy time series. 

This paper is organized as follows: Section III] is devoted 
to the noise-free stochastic process, the definition of its joint 
probability density and expressions for the moments of this 



density in terms of a Taylor-Ito expansion. Section III pro- 
vides the properties of the measurement noise under consid- 
eration and in section |IV] expressions for the moments of a 
noisy process will be derived. After looking at the benefits 
of equidistantly sampled experimental time series in section 

V] the previously derived expressions will be used in section 

VI] to extract the parameters of the measurement noise and in 
section VII to extract polynomial approximations for drift and 



diffusion functions. Finally in section VIII the results will be 
applied to some synthetic time series. The used properties of 
the Gauss function and further computational details are given 
in appendices [A] and [B] 



II. STOCHASTIC PROCESS 

Let X(i) be a stochastic process in N dimensions that can 
be described by a time-independent Ito -Langevin equation 



dXi{t) = D(i)(X)df+ JD(2)(X)dW(i), (1) 



where D^^^ and T>'-^^ are the Kramers-Moyal coefficients of 
the corresponding Fokker-Planck equation and dW denotes 
a vector of increments of independent Wiener processes with 
< dWidWj >— Sijdt. The notation vD(^ is used to denote 
a matrix g with g • g* = D^^^ IJSTl . 

Let the one- and two-point probability density functions of 



X be denoted by 



p(x) := p(x,<) 
p(x,x',t) := p{:>i,t;x,t+T) 

= p(x)p(x',i+T|x,t) 



(2a) 
(2b) 
(2c) 



and let the conditioned moments of p{x, x',t) be defined 
as follows (the notation dx is used to denote the product 
dxi . . . dxN whereas dx denotes a vector of differentials dxi). 



m(")(x) = I p{x,x',T)dx' 
ml'{x,T) = {xi-Xi)p{x,x.',T)dx' 

J-x.' 

"^1j(x,t) = {x'i-x,){x'j~Xj)p{x,yi',T)dx' (3c) 



(3a) 
(3b) 



These moments are observable quantities. For a given time 
series they can be estimated by binning or other density- 
estimation techniques. Using Eq. (|2c]l allows to express the 
moments m'^'^) in terms of moments h^*^' of the conditional 
increments of X 



m(")(x) = p(x)-l 

"^r'(x,T) = p(x) •/lf^(x,T) 

m'^\x,r) = p(x)./z(f(x,r), 



(4a) 
(4b) 
(4c) 



The explicit terms up to second order are given below (using 
index notation and summation convention). A detailed de- 
scription of the Taylor-Ito expansion and its moments can be 
found in Bm. 
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D^d^D. 



(1) 



(7a) 



h)V = tD. 



(2) 
i(2); 



T 



i(l)n(i) 



2Dr>D\'>+D)VdkD 



(2). 



l(l) 



ik 



^D);'dkD 



(1) 



-dI}'^o,dS^ 



■^n(2);) f, n(2) 



+ 0{t') 



(7b) 



Inserting the series representations into Eq. (Hll yields a rela- 
tion between the observable moments m'^'^' and the unknown 
functions D^^^^ and D'^2) por small values of t this allows the 
direct estimation of the Kramers -Moyal coefficients. 



i^f^(x) 



<'(x) 



1 mf'{x,T) 

T to('^)(x) 

lm|f (x,r) 



^(0)(x) 



0{t) 
0{t) 



(8a) 
(8b) 



III. MEASUREMENT NOISE 



with 



h^\x,T) := <[X,{t+T)-X,{t)] 



X(t)=x 



> 



/*!f(x,r) 



:= <[X,{t+T)-X,{t)] 



x[X,{t+T)-X,{t)] 



x(t). 



> 



(5a) 



(5b) 



A Taylor-Ito expansion of Eq. ([T]i provides expressions for 
these expectation values. Provided that D^^^ and D(2) ^re 
smooth functions in x, it is possible to represent h'^^ and h(2) 
as power series in r. The lowest order terms in these series are 
linear in r. The series-coefficients are given by sums of prod- 
ucts of the Kramers-Moyal coefficients and their derivatives 
and are thus generally functions of x. 



h^\x,r) = ^c^^'Hx)r'^ 



fe=i 



^i?(-'-) = E^&'^w-' 



(6a) 
(6b) 



fc=i 



The measurement noise under consideration is denoted by 
Y(i) and described by an Omstein-Uhlenbeck process in N 
dimensions. Such a process is characterized by linear drift- 
and constant diffusion functions and its statistical properties 
can be derived analytically (see e.g. |31|). The temporal evo- 
lution of Y is described by Eq. (|9]l. Here the eigenvalues of 
matrix A are required to have positive real part and matrix B 
is assumed to be symmetric and positive semi-definite. The 
notation \/B is used to denote a matrix g with g • g* = B 
and the elements of dW denote the increments of indepen- 
dent Wiener processes with < dWidWj >= Sijdt. 



(9) 



dY{t) = -A-Ydt + VBdW{t) 



While A and B are appropriate to describe the temporal evo- 
lution of Y, the 'macroscopic' properties of the noise are 
more conveniently described in terms of the covariance ma- 
trix V and the matrix of (exponentially decaying) correlation 
functions M(r). Looking at the auto-covariance of Y one 
finds 



<Y{t+T)Y\t)> = M(t)V 



(10) 



with 



M(t) 



„-At 



V = 



^^'Be-^''ds. 



(11) 
(12) 



Furthermore Y is found to be Gaussian distributed. If 
G(V, x) is used to denote a normalized Gauss function in x 



with covariance V (see Eq. ( A4 1), then the one- and two-point 



probability density functions py can be written as 



py(x) 
Pi'(x,x',t) 

with the shortcut 



G(V,x) (13a) 

G(V,x)G(C(T),x'-M(r)x) (13b) 



C(t) := V-M(r)VM*(T). 



(14) 



A note on the eigenvalues of A and M: If A^ denotes the 
eigenvalues of A, then the eigenvalues of M(t) are given by 
g-AiT gy introducing the relaxation times (or characteristic 
time scales) Ti as Ti :— 1/Xi, the eigenvalues of M can be 
written as e""^/^'. In the numerical examples given later, the 
measurement noise will be characterized by such relaxation 
times Tj instead of by eigenvalues of A. 



IV. NOISY STOCHASTIC PROCESS 

Let X*(t) = X(t) + Y(t) be the sum of a stochastic sig- 
nal X(t) and measurement noise Y(i) as introduced in sec- 
tions III] and III respectively. Because X and Y are indepen- 
dent stochastic variables, the probability density functions of 
their sum X* is given by the convolution of the individual 
density functions. 



p* (x, x', t) = Py (x, x', r) * p(x, x', r) 

= / / /9Y(x-z,x'-z',r) 
X p(z, z', r) dzdz' 



(15) 



Instead of the conditioned moments m^'^) only their noisy 
counterparts m*''''^ can be determined. 



m*(")(x) 



m^ '(x,r) 



m^ (x, r) 



/ /9*(x,x',r)dx' 

j {x'i~Xi)p*{x,x',T)dx' 



(^i Xi)(Xj Xj) 

xp*(x, x',r)dx'. 



(16a) 
(16b) 

(16c) 



integration with respect to x'can be performed within the con- 
volution integral. Using the definition of the moments m^'^) 
and taking advantage of the properties of the Gauss function 
then finally leads to the following equations (see appendix [B|. 
Function arguments are omitted for notational simplicity. 



*(i) 

*(2) 



,(0) 



Py * rn 



(17a) 
(17b) 



+Qii'din 



*(i) 



Qjfd,,m\ 



(1) 



(17c) 



Here Py(x) — G(V,x) is the density function of the mea- 
surement noise Y. The terms h'^''), as introduced in section|n| 
denote the moments of the conditional increments of X. The 
quantity Q, finally, has been introduced as an abbreviation and 
is defined as 



Q(r) 



(Id - M(t))V. 



(18) 



Equation ( 17 1 allows to express the observable moments 
m*'^'^' in terms of the unknowns rrS^\ h.^ \ M and V. How- 



ever it is possible to use Eq. (17bi to extract the parameters 



of the measurement noise without the need for a simultaneous 
determination of h^''' and ttt,^^^ This will be done in section 
VI Next, however, an assumption on the given time series will 
be made. 



V. EXPERIMENTAL TIME SERIES 

It will be assumed, that the values of a given time series are 
taken at equidistant points in time with a basic time increment 
of Ai. This is often assumed tacitly but shall be stated here 
explicitly because it will be used in the following. 



X*(iAt), 



1, 



(19) 



Increments of X* can thus be calculated for all time incre- 
ments T which are integral multiples of At. The experimental 
two-point probability density p*(x, x', r) for those values of 
T can then be written as a sum of Dirac -distributions. 



p*(x,x',fcAt) = 



1 

X rC 

5(x' 



E 

X 



<5(x-x:) 



i+k) 



(20) 



Weighted integrals of the 'true' density p* can easily be esti- 
mated by weighted integrals of p* now. Given a weight func- 
tion /(x, x') and denoting the estimate by / one finds 



Inserting the definitions of p* and py (Eqs. (15 1 and (13bi 
respectively) and interchanging the order of integration, the 



/ = 



/(x, x')p* (x, x', kM) dx dx' 



k ^ 



/(X*,X*_,_j,). 



(21) 



integral power of M(At) then. This is due to the fact that 
M(r) (according to section III i is a matrix exponential. 



Weighted integrals of p* can therefore directly be estimated 
from the given time series. There is no need to use binning 
to estimate the density p* first |29|. Weighted integrals of the 
moments m*^*^) can be treated the same way by expressing 



them as weighted integrals of p* using Eq. ( 16 1 



Mo 
^ M(fcAi) 

So finally one gets 



-AAt 



:= M(At) = e 



(27a) 
(27b) 



VI. EXTRACTING MEASUREMENT NOISE 
PARAMETERS 



Multiplying Eq. ( 17b i by Xj and subsequently applying an 



integration with respect to x leads to 



*(i) 



dx = f[pY*ih\^^m^°'>)]xjdx 

+Q^^' f {d^,m<°^)x,dx. (22) 



The left hand side of this equation can directly be estimated 
from a given time series and will be denoted by Z. 






(23) 



Applying integration by parts allows the evaluation of the sec- 
ond integral on the right hand side 



Qie / {d,:m*^''^)xjdx 



(24) 



The remaining integral in Eq. (22 1 only depends on t because 
of the function h'^^)(x, r). Using Eq. (6ai therefore allows 
to express the integral as a power series in r with unknown 
coefficients P^ ■ Truncating this series to some order i/^ax 
yields an approximation of the integral by a polynomial in t. 



[py*{h''Pm(^))]x,dx ^ Y.W^ 



i^)^v 



(25) 



Putting this all together (and additionally replacing the abbre- 
viation Q by its definition) so far yields 



Z(fcAi) = ^p('')(fcAt)'^-(ld-M^)V. (28) 



Evaluating Z for a sufficient number of increments, kAt, 
yields a system of equations that can (iteratively) be solved 
for the unknowns P'"^', V and Mq in a least square sense. 
Subsequently the relaxation times of the measurement noise, 
Ti, can be calculated from the eigenvalues of Mq. 
For such a fit to succeed, two conditions must be met. 

• Firstly, the largest increment kj„.^j^At should be small 
compared to the characteristic time scale of the under- 
lying stochastic process. This will allow, to chose a low 
polynomial order j/max (the smaller r the better h'^^ can 
be approximated by low order polynomials). 

• Secondly, the relaxation times T^ should be small com- 
pared to fcmaxAt. This will allow, to distinguish the ex- 
ponential functions in M(r) from a low order polyno- 
mial. 

The proposed method is therefore limited to measurement 
noise with relaxation times T, considerably smaller than the 
time scale of the underlying stochastic process. 



VII. EXTRACTING DRIFT- AND DIFFUSION FUNCTIONS 



In the following it will be assumed, that the noise parame- 



ters have already been estimated according to section VI The 
parameters and derived quantities like Q(t) will therefore be 
treated as known quantities. 



Multiplying Eqs. ( 17b i and ( 17c i by some weight function 
vl'(x) and subsequently applying an integration with respect 
to X yields 



Ihs, 



* 



Ihs,, = / * 



PY * (/ir^mt")) 



dx 
dx 



(29a) 
(29b) 



Z(t) = ^P^'^V''- (Id-M(T))V. 



(26) 



where Ihs, and Ihs,, are abbreviations for the left hand sides 



Assuming that the time series is sampled with a basic time 
increment Ai (as stated in section fv\, the value of Z(r) can 
directly be estimated for all increments r being integral mul- 
tiples of At. The corresponding value of M is given by an 



Ihs, 



* 



Ihs,, := / * 






dx (30a) 



(Q 



ij 1^ ^ji 



-QH'Qjrd,>dj.)m<°^ 



dx. (30b) 



Applying integration by parts allows to express the integrals 
in Eq. ( 30 1 as sums of weighted integrals of m*^*^) . For exam- 
ple one finds / *to*'^^'' + Qu^ J {di'^)m*^"'> for the left hand 
side of Eq. (|29a|. This expression can directly be estimated 



from the given time series because Q and ^ (and thus also the 
derivatives of ^) are known. The same holds for the left hand 
side of Eq. ( |29b| l. Both left hand sides can therefore directly 
be estimated for a given choice of t and ^P. 

The corresponding right hand sides, however, refer to the 
unknown function rn(°\ It is possible to overcome this prob- 
lem if the drift- and diffusion functions are approximated by 
polynomials in x. 



i(i) - 



D^? = a. 



^2a ' 



DT' = aT' + ay^Xa + a)^;pXaXfi 






, (2) , (2) 



^ija-^a 



\jal3-' 



(31a) 
(31b) 



The coefficients in the power series representation of the con- 
ditional moments h^'^'' then also become polynomials in x. 



with the linear differential operators 

Lap = Vaa'Vpp'da'dp' —Vap 



(35) 



Applying integration by parts then allows to express Eq. (33 1 
in terms of weighted integrals of m*(°\ which can directly be 
estimated from the given time series. One finds 



F = *m*(°)dx 
Fa = f {*a;„m*(°) - [L„*]m*(°)} dx 



(36) 



Expressing the right hand sides of Eq. (29i in terms of 
F,Fa,--. one finally obtains the following equations. 



.(1) - 



.(2) 



(1) , (1) , (1) 

'b?''+b^jxa + b^jpxaxp- 



,(2) 



l(2) 



,(2) 



,(2) 



~ij + ^ija^a + (^ijal3^aXp 



eixa 



(.(2) 



'ij ' ^ija-^a -r bj^j'^i^XaXjS 



(32a) 



(32b) 



Ihs,; 



+ T 
+ . 



lhs„ = 



+ r 



(2)p , „(2) 






a 



uC^) p , l(2) p , . (2) p , 



^'^ F , 

■'^ 1 
ijaP 



(37a) 



(37b) 



For the sake of simplicity the abbreviations b^'") have been 
introduced here. However all the coefficients in Eq. ([32] i can 
of cause be expressed in terms of the coefficients a^*^-'. Us- 



ing Eq. ( 32 1 the problem of expressing the right hand sides 



of Eq. ( 29 1 reduces to the problem of expressing terms of the 
form 



Evaluating the left hand sides and the quantities F,Fa,-.- for 
a sufficient number of increments r and weight functions ^ 
yields a system of equations that can be solved for the un- 
known polynomial coefficients in a least square sense. There 
are different approaches to deal with the higher order terms in 
r now. 



Fr, 



:= vp 



PY * {Xai ■■■Xaf.m^"'') 



Because py = G(V, x) is a Gauss function, the convolution 
within the square brackets can be expressed in terms of deriva- 



tives of TO*(°', Xam*^°\ XaXpm*'^'^^ 

One finds 



(see appendix 



A5i 



pY * [Tin 

Py * [xafin 

Py * [xaXpm 



(0)] ^ 
(0)1 _ 






T *(0) 

XaXfjin*^'^' + La[xpm*^^'] 



• The most simple approach is, to completely ignore the 
dx. (33) higher order terms. This will lead to a linear fit in r. 

Furthermore the resulting set of equations will be linear 
in the unknown coefficients a'^'''. 

• A more elaborate approach is, to perform a polynomial 
fit in T. Coefficients beyond some order will be ignored. 
If the remaining coefficients are treated as additional 
unknowns, then the resulting set of equations will stay 
linear However this way some available information is 
ignored, because b*^*^^ and higher order coefficients can 
in fact be expressed in terms of the coefficients a''^^ 

• Finally, a polynomial fit in r can be performed, where 
the higher order coefficients are expressed in terms of 

(34) the coefficients a'*^'. This will improve the accuracy 



of the estimate, because of the smaller number of un- 
knowns. The resulting set of equations, however, will 
now be nonlinear and needs to be solved iteratively. 

The choice of the weight functions $ has been left open up 
to now. Obviously ^P needs to admit the various integrations 
by parts that have been applied. For these integrations it also 
has been tacitly assumed that the involved boundary values at 
|x| — ?► oo are vanishing. This imposes additional restrictions 
on ^. 

The set of weight functions used in the numerical examples, 
given below, consisted of a number of Gauss functions cen- 
tered at different points in space. There may be better choices, 
but the problem of finding the 'best' set of functions will not 
be addressed here. Gauss functions are smooth and real val- 
ued and have a local support. But maybe it would be better 
to choose, for example, some complex valued functions like 
exp(ztL'*x) which have a local support in Fourier space only. 
Also piecewise polynomial functions like the B-spline base 
functions may be an alternative. 



VIII. APPLICATION TO NUMERICAL DATA 

In order to check the accuracy of the proposed method, a 
test case in two dimensions has been investigated. A stochas- 
tic process, X(t), as introduced in Sec. Ill] has been specified 
by the following choice for the drift- and diffusion functions 
(x and y denote the components of the vector x). 



The measurement noise, Y(t), as introduced in Sec. Ill is de 



D(i)(x) 
D(2)(x) 



X — xy 
x^ -y 



0.5 

0.5(1 + a;2) 



(38a) 
(38b) 



By numerical integration synthetic time series of the process 
can be generated. All series used in the following will consist 
of 10^ points, sampled at time increments At = 0.005. The 
deterministic part of the process dynamic and the experimen- 
tal density distribution of X is visualized in Fig.[T| 




scribed by an Ornstein-Uhlenbeck process in two dimensions. 
The noise is characterized by eigen-directions, e^M. and cor- 
responding relaxation times, Ti, of its matrix M and by the 
principal directions, e^v, and the corresponding standard de- 
viations, CTi, of its covariance matrix V. The following values 
have been chosen (relaxation times are given in units of At). 



(ei,e2)M 
(ei,e2)v 



1 1 

2 J ' 

1 -1 
1 1 



0.25 
0.5 



(39a) 
(39b) 



The deterministic part of the dynamic of the measurement 
noise and the experimental density distribution of Y is visu- 
alized in Fig. [2] 



b) 




12 



FIG. 2: Deterministic dynamic and density distribution of the 2D 
measurement noise Y(f). The trajectories in phase space (a) are 
generated by the deterministic part of the process dynamic (x = 
—X + y/3, y — ~y/3). There exists a single, attractive, fixed point 
at the origin. The contour lines (b) of the probabihty density function 
have been computed from a time series of Y. 

Adding the time series of X and Y yields a series of 'noisy' 
values X*(i) = X(i)+Y(t). This will be called a noisy time 
series in the following. Excerpts of X(i) and X*(i) as well 
as the experimental density distribution of X* are shown in 
Fig. [3] 



a) 



l^py^fMfi^^ 





-3 -2 -1 



FIG. 1: Deterministic dynamic and density distribution of the 2D 
process X(f ). The trajectories in phase space (a) are generated by the 
deterministic part of the process dynamic (x = x — xy, y = x^ — y). 
There exist three fixed points: a saddle at the origin and two stable 
foci at ix = ±l,j/ = 1). The contour lines (b) of the probability 
density function have been computed from a time series of X. 



FIG. 3: Excerpt of a time series of X(i) (a) and of a corresponding 
noisy time series of X*(i) (b). The contour lines (c) of the prob- 
ability density function have been computed from of a noisy time 
series. 

Now the extraction of the measurement noise parameters, 
as described in Sec. VI has been tested. For a sample of 



1000 independent realizations of noisy time series the matri- 
ces M(Ai) and V have been estimated. For each estimate a 
number of scalar quantities has been calculated. For a char- 
acterization of the deterministic part of the noise-dynamic the 
angles, ai, spanned by the eigendirections of M and the x- 
axis, and the relaxation times, Ti, determined by the eigen- 
values of M(Ai), have been used. Their true values are given 
by 



ai ^ 0°, a2 ~ 63.43° 
Ti/At = 1, Ts/At = 3 



(40a) 
(40b) 



The covariance matrix V is symmetric and can thus be char- 
acterized by three scalars: the angle /3, spanned by the first 
principal direction of V and the x-axis, and the standard de- 
viations (Ti in direction of the principal axes. The true values 
are given by 



/3 

CTl 



45° 
0.25, 



0-2 



0.5 



(41a) 
(41b) 



Parameter fitting has been performed with a maximum in- 
crement Tmax = 50Ai and a maximum polynomial order of 
t'max = 3. The resulting distributions of the estimates are 
shown in Figs. HI and l5] It turns out, that the sample standard 
deviations of the angular quantities are given by some tenth of 
a degree. Relaxation times and noise strengthes are estimated 
with relative errors well below one percent. 




99 0.995 1 1 005 1.01 



FIG. 4: Histograms of the observed distributions of the estimated 
parameters characterizing matrix M of the 2D measurement noise. 
Angles Qi of the eigen-directions (a,b) and corresponding relaxation 
times Ti in units of Ai (c,d). The standard deviation of the respective 
distribution is given by an annotation. The true parameter values are 
indicated by solid vertical lines. 

For the estimation of the drift- and diffusion functions a com- 
plete quadratic ansatz has been made for each component of 

D^^) and D*^-^'. Because D^^^ is symmetric, this leads to a 



a) 




0.498 0.499 



0.501 0.502 



FIG. 5: Histograms of the observed distributions of the estimated 
parameters characterizing the covariance matrix V of the 2D mea- 
surement noise. Angle /3 of the first principal direction (a) and the 
principal standard deviations ai (b,c). The standard deviation of the 
respective distribution is given by an annotation. The true parameter 
values are indicated by solid vertical lines. 



total of 30 coefficients. As maximum time increment for the 
fitting procedure a value of Tn,ax = 15Ai has been chosen. 
The set of weight functions 'J consisted of 16 Gaussian func- 
tions centered at the nodes of a rectangular 4x4 grid cover- 
ing the ±2(7 range of the experimental density distribution of 
X* . The standard deviations of the weight functions itself was 
chosen as twice the distance between neighbouring nodes. 

For this setup the coefficients have been estimated now for 
a sample of ten independent realizations of the noisy time se- 
ries. Using a linear fit in r leads to the results shown in Figs.l6] 
and|7] 

It can be seen that some of the estimates, especially for the co- 
efficients of the diffusion functions, are significantly biased. 
Looking, for example, at coefficient 02 of diffusion function 
D[^' one finds a value of —0.1663 ± 0.0097 which signifi- 
cantly differs from the true value of zero. Much better results 
are obtained by performing a quadratic fit in t. To do so, the 
most simple approach has been chosen: For each parameter 
Oq an additional parameter b^ (see Eq. ( 37 1) is introduced. 



The only purpose of this parameters is, to absorb some of the 

quadratic terms in r. Because the number of unknowns is 

doubled this way, this will also lead to higher fluctuations of 

the estimates. However, performing such a quadratic fit also 

greatly reduces their biasing, as can be seen in Figs. 18] and 19] 

(2) 
For the above mentioned coefficient 02 of D\^^ 

obtaines a value of 0.0025 ± 0.0352. 



.e.g.. 



The extraction of noise and process parameters from a noisy 
time series seems to work for the given 2D test case. To check 
if this also holds for higher dimensions, the test case has been 
extended to four dimensions. Process and measurement noise 
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FIG. 6: Parameter estimates for the drift functions D[ (a) and Dj 
(b) of tlie 2D process, obtained by a linear fit in t. For each polyno- 
mial coefficient the calculated estimates are given by grey bars with 
an additional white bar to the right, which shows the true value. 



FIG. 8: Parameter estimates for the drift functions Dl (a) and Dj 
(b) of the 2D process, obtained by a quadratic fit in r. For each 
polynomial coefficient the calculated estimates are given by grey bars 
with an additional white bar to the right, which shows the true value. 
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FIG. 7: Parameter estimates for the diffusion functions D\i (a), FIG. 9: Parameter estimates for the diffusion functions D\^ (a), 

(2) (2) 

DJ2 (b) and D22 (c) of the 2D process, obtained by a quadratic 
fit in r. For each polynomial coefficient the calculated estimates are 
given by grey bars with an additional white bar to the right, which 
shows the true value. 



(2) (2) 

Dj2 (b) and D22 (c) of the 2D process, obtained by a linear fit in 

r. For each polynomial coefficient the calculated estimates are given 

by grey bars with an additional white bar to the right, which shows 

the true value. 
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now are defined in x, y, z, w space. The process is defined by 



D(i)(x) 



D(2)(x) 



fx — xy\ 



—z 
—w 



n 



ii^ 
" 2 








2 



Vo 

and the measurement noise by 






2 ' 



(61,62,63,84) 



M 



(ei,e2,e3,e4)v = 



'1 1 1 ON 

2 

2 

^0 ly 

a -1 ON 

110 

10 

vO 1/ 



(42a) 



(42b) 



(43a) 



(43b) 



weight functions $ consisted of 81 Gaussian functions cen- 
tered at the nodes of a rectangular 3x3x3x3 grid covering 
the ±2(T range of the experimental density distribution of X* . 
The standard deviations of the weight functions itself was cho- 
sen as twice the distance between neighbouring nodes. Using 
the same type of quadratic fit in t as in the 2D case yields the 



coefficient estimates shown in Fig. 10 



D1J 

D1_2 

D1_3- 

D1_4 

D2_11 

D2_12 
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D2_22 

D2_23 

D2_24 

D2_33 

D2_34 

D2_44 



a3 a4 a11 a12 a13 at 4 a22 a23 a24 a33 a34 a44 
polynomial coefficient 




(44) 



Estimating the relaxation times and the noise strengthes for a 
noisy time series (10^ points. At = 0.005, r^ax = 50, I'max — 
3) yields the following results. 
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FIG. 10: Parameter estimates for the polynomial coefficients of 
the drift- and diffusion functions of the 4D process, obtained by a 
quadratic fit in r. For each coefficient the calculated estimate is 
given by a black bar. White bars show the respective true values. 
Black horizontal lines show the abscissa of the coordinate systems 
of the respective Drift- and Diffusion functions denoted to the left. 
The results for constant-, linear- and quadratic coefficients as well 
as the results for drift- and diffusion functions are seperated by grey 
lines. 



IX. CONCLUSIONS 



The accuracy of the estimated matrices, M and V, can be 
expressed in terms of the relative errors cm and ev, defined 
as M^^M — Id and V^^V — Id respectively. One finds 



Cm 



ev 



-0.6 -5.2 -5.2 +2A\ 

-0.2 -F0.8 -f2.2 -1.4 

-0.5 +0.8 -3.0 +0.8 

^-1.0 -1.5 -F2.0 -1.5y 

^-hO.7 -0.6 -fO.7 -1.7\ 

-2.1 -h3.1 -1-1.8 -2.5 

-0.8 +3.3 -5.0 +1.9 

-0.5 -3.8 +1.9 -1.4/ 



X 10" 



X 10 



-3 



(46a) 



(46b) 



For the estimation of the drift- and diffusion functions a com- 
plete quadratic ansatz has been made for each component of 
D'^-* and D'^'^-'. In four dimensions this leads to a total of 210 
coefficients. As maximum time increment for the fitting pro- 



cedure a value of Trr 



15 At has been chosen. The set of 



A procedure has been described for the analysis of stochas- 
tic time series in N dimensions in the presence of strong mea- 
surement noise. The algorithm is able to cope with exponen- 
tially correlated noise and accurately extracts strength and cor- 
relation time of the measurement noise as well as the parame- 
ters defining the drift- and diffusion functions of the underly- 
ing stochastic process. This has been shown by the analysis of 
synthetically generated time series in two and in four dimen- 
sions. 

The ability to deal with exponentially correlated measure- 
ment noise in more than one dimension has not been given by 
the approaches available up to now. 

Because of the use of weight functions there is no need to 
perform any density binning. All required quantities can be 
obtained from weighted sums of the values of the time series. 
This avoids the aliasing errors caused by finite bin sizes. 

All calculation have been performed on a standard desktop 
PC. The analysis of a signal took about five minutes (2D case) 
respectively fifty minutes (4D case). 
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In the current implementation only a simplified quadratic 
fit in the increments r can be performed. Implementing a full 
polynomial fit, as mentioned in Sec. VII should allow to ex- 



tend the range of time increments that can be used for the 
analysis and thus should increase the accuracy of the results. 
This has to be done in the future. 

Another point to be improved is the restriction on polyno- 
mial approximations of the drift- and diffusion functions. An 
approximation by spline-based functions would be much more 
flexible. When using such a parametrization, however, it will 
no longer be possible to accurately express the convolutions 
in Eq. (29 1 in terms of observable quantities. It will become 
neccessary to also introduce a parametrization for the density 
77i(o) which significantly complicates the calculations and also 
introduces additional parameters to be estimated. 

Also a future task is the application to some real world data. 
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Above definitions imply the following properties. 



•^[^/(^)] = *^./('^) 



dx 



^[/(x)*.g(x)] = 



_d_ 
dxi 



[/(x)*g(x)] 



f{uj)g{i^) 
, d 



dxj 
= /(x) 



/(x)] * .g(x) 



dx^ 



5(x)] 



(A3a) 

(A3b) 
(A3c) 

(A3d) 



3. Derivatives and monomial products of Gauss functions 

Let G(C, x) denote a normalized Gauss function with co- 
variance matrix C and function argument x G M^. 



G(C,x) := 



ze 2 



ix'C-^x 



Appendix A: Gauss functions 



V(27r)^|det(C 
This function is a eigenfunction of the Fourier transform 



(A4) 



1. Index-vectors and monomials 

For the sake of a compact syntax, multiple indices will 
frequently be combined into an index-vector. For example 



G{C,uj) = e-5 



— TjUJ i^UJ 



A, 



will be written as A; . To denote the length of such a 



vector j, the function ^(j) will be used. 

As a further abbreviation the symbol A4 is introduced 
for monomials of the components of a vector. A mono- 
mial Xjj . . . Xj^^ will be written as Mji...j„ (x) or simply as 
Mj (x). Monomials of the nabla vector will be used, to denote 
multiple partial differentiation more compactly by Alj(V). 



2. Fourier transform and convolution 

Let Fourier transform and convolution of functions M.^ — > 
C be defined as below. For notational simplicity the 'hat' 
syntax will be used to denote the Fourier transform of sin- 
gle functions. For more complex expressions the functional 
form J^(. . .) will usually be the better choice. 



^[/(x)](u;) = /(a;) := / e"*'^ V(x) dx (Al) 
/(x) * .9(x) := / /(x').g(x - x') dx' 

Jyi' 

= j /(x-x')ff(x')dx' 



(A2) 



= ^(27r)^|det(C-i)|G(C-\u;) (A5) 

It can be shown by mathematical induction, that the deriva- 
tives of G all have the form 



7Wj(V,)G(C,x) = Pj(C,x)G(C,x), 



(A6) 



where P is a polynomial of order £(j) in x. Mathematical in- 
duction also shows, that P only contains monomials in x of ei- 
ther even or odd order. The coefficients of P can be expressed 
in terms of the elements of C^^, but no attempt will be made 
here to give an explicit formula, because the expressions for 
any finite order polynomial can be derived iteratively. Up to 
order three the derivatives of G are given by (using summation 
convention) 



dxj 


— Gj^ Xa 


G 






(A7a) 


^" G - 
dxjdxk 


^ja^kP ^aXp — 6^- J, 


G 


(A7b) 


^' G - 
dxjdxkdxi 


— Gj^Cf.pGi^ XaXpXj 

HCi^Cr^+Gj-^G^'^ 






+Cklc-, 


I) -. 


G. 




(A7c) 
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Terms of the form Ai{x)G will be called monomial products 
of G in the following. Such products can be expressed in 
terms of derivatives of G. Applying a Fourier transform to 



Eq. ( A6 1 and using of Eqs. ( A3a I, ( A3b i and ( A5 1 first gives 



7Wj(ia;)G(C-\a;) = Pj(C, iV<^)G(C-\a;). (A8) 

Substituting iu! by x and C by — C^^ (thus u)*'Cuj by 
x*C^^x) then finally yields 



Xj(x)G(C,x) - Pj(-C-\-V,)G(C,x). (A9) 
Up to order three the monomial products of G therefore read 



5. Gauss functions in convolutions 

Convolutions of the form [A^j(x)G(C,x)] */(x) can be 
expressed in terms of derivatives of the convolution G*f. This 
can be derived straightforwardly by first expressing A^(x)G 



by derivatives of G and then applying Eq. ( A3d i. One finds 



[Xj(x)G(C,x)]*/(x) =. Pj(-C-i-V,) 

x[G(C,x)*/(x)]. (A12) 
It is also possible to express convolutions of the form 
G(C, x) * [A^j(x)/(x)] by derivatives of monomial products 
of G* /. This can be derived in Fourier space. So let F denote 
the Fourier transform of the expression under consideration. 



XjG 



iAj o Jj ji^ v_T 



XjXkXiG 



Cja 



d 

dXr. 



G 



(^ja^-^k/S 



d^ 



dxadxp 

(-'jaGkpGlj ■ 



- Cjk 

g3 



(AlOa) 
G (AlOb) 



dxadxpdxj 

'[(-yjkGla + GjlUka 

d 



~\~Gkl(-'ja] 



dXa 



G. 



(A 10c) 



F := ^{g(C,x)*[Xj(x)/(x)]} 
= G(C,a;)Xj(zV^)/(a;) 



(A13) 



Using the identity 



1 



G(C,a;) 
leads to 



= A/(27r)^|det(C-i)| G(C-^ia;) (A14) 



4. Moments of Gauss functions 



Integrating Eq. ( A9 1 with respect to x yields expressions for 
the moments of G. Because integrals of derivatives of G are 
vanishing, the moment is determined by the constant part of 
the polynomial Pj . This coefficient will be non-zero only for 
even moments. The odd moments of G all evaluate to zero (as 
can also bee seen from symmetry considerations). The first 
non-vanishing moments are given by 



F = G(C,u;)Xj(.V.)[g^^/M] 

G{C,uj) 



= G(C,u;)^(2^)^|det(C-l)|Xj(^V„) 
x{G{C-\iuj)[G{C,uj)f{cj)]}. 



(A15) 



Now the product rule of differentiation is applied to the term 
in the curly brackets. Using index- vectors the product rule can 
be written as 



/ Gdx 


= 1 


(Alia) 


Jx 






/ XjXkGdx 


= Cjk 


(Allb) 


Jx 






/ XjXkXlXmGdx 


— CjkClm + CjlCkm 




Jk. 


+ CjrnCkl- 


(Allc) 





A^j(V)[/(x)5(x)] 



J2 [A^j'(V)/(x)] 

(J'J")6-P(J) 

x[Xj„(V)5(x)]. (A16) 



Here V{j) denotes the set of all 2^^^'> pairs (j', j") that can be 
obtained by distributing the components of j on two vectors j' 
and j". Applying the product rule yields 



F = G(C,a;)y^(2^)^|det(C-i)| ^ {Xj-(zV„) G(C-i, zu;)}{Xj-,(jV^)[G(C,u;)/M] }■ (A17) 



(j'J")e7'(j) 
Temporarily substituting z = —iu) in the first bracket and using G(», — x) = G(», x) gives 



My{iS7u,)G{C-\iLj) = My{W,)G{C-\z) = Py{C-\z) G{C-\z) = Py{C-\-iu:) G{C-\iuj>). (A18) 
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Now G{C \ io;) can be written in front of the sum. Using Eq. ( A14 1 some factors cancel out and it remains 



(J'J")6T'(j) 

Switching back to real space finally gives the desired relation 



(A19) 



G(C,x)*[Xj(x)/(x)] - Y. ^j'(C-^-V,){Xj.(x)[G(C,x)*/(x)]}. 



(A20) 



(j'J")eT'(j) 



(k) 



Appendix B: Conditioned moments m 

The somewhat lengthy calculations leading to Eq. ( [I7| ) are given below. The function argument of M(T)and C(t) will be 
omitted for notational simplicity. Partial derivation with respect to Xj will be denoted by 9, and Einsteins summation convention 



will be used. Starting with Eq. (16 1, inserting Eqs. (13b i and ( 15 i and interchanging the order of integration gives 



TO*(°)(x) = f G{V,x-z) f p{z,z',t) f G{C,x'^z'-M-{x-z))dx'dz'dz (Bla) 

Jz J z' Jx.' 

m*^^^(x,r) = f G{V,x-z) f p{z,z',t) f {x'^~x^)G{C,x'-z'-M■{x-z))dx'dz'dz (Bib) 

Jz J z' Jx.' 

'(x,r) = I G{Y,x-z) I p{z,z',t) I {x'i-x{){x)-Xj)G{C,x'-z'-m-{x-z))dx' dz' dz (Blc) 

Jz J z' Jx' 



*(2)/ 



The integrals with respect to x' can be expressed in terms of the moments of the involved Gauss function (see Sec. (A4l). Using 
the definition of m*^''^ then allows to express m**^") by a convolution. 



(0)( 



*^")(x) = G(V,x)*to(")(x 



The other moments so far read 



*(i) 



(x,r) = / G(V,x-z) / p{z,z',T)[zl-Xi+Mii^{xi^-Zi,))] dz' dz 

J Z J z' 

m*f\x.T) - I G{Y.x-z) j p{z,z'.T)[C,, + {z[~x, + Mu'{x,,-z,,)) 



X [z'a — X j + AI j ji {x j' — Zj'))] dz' dz. 



(B2) 
(B3) 



(B4a) 



(B4b) 



Sorting the terms in rectangular brackets by powers of the components of z' — z and using the definitions of the moments m'^'^-' , 
allows the intergals with respect to z'to be expressed by the moments m'^'"'^. 



m*^'^(x,r) = /g(V,x-z) m^^\z,T)-{5u'-^■Ue){x^>-z,,)m^''\z) dz 

J Z 

mf\x,T) = /'g(V,x-z) \m^;^\z,T) + C.,m<-°\z)-{5,y-Mj,,){xy^zy)m^p{z,T) 

- {5ii> - Mii> ) {xi> - Zi, ) my' (z, t) + ((5ij/ - M^i, ) (a;^/ - z^/ ) {Sjy - Mjy ) {xy -zy)m^^\z) 



(B5a) 



dz. (B5b) 



Now the relation / [xi — Zi)f{x — z)g(z) dz — [xif{x}] * g{x.) can be used to express the noisy moments as convolutions. 
Function arguments can now be omitted without confusion (G refers to G(V, x)). 



*(i) 

i 
*(2) 



,(1) 



(0) 



V 



G*m\ ' - ((5ij' - Afjj.)[a;j'G] * m 

G * m^f + G,, G * m(o) - {Su'-Mu')\x,,G] * m^ _ {6j^,-Mjj,)[x,>G] * 



(B6a) 



,(1) 
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+{Sii>-Mii'){Sjj>-Mjj>)[x^'Xj>G] * m 



(0) 



(B6b) 



Because G is a Gauss function, the monomial products XiG and xiXjG can be expressed by derivatives. Inserting the definition 



of C (Eq. ( 14 1) and resorting terms subsequently leads to 



m*« = G^mf^ + {5,^,~Nh,,)V,,k[dkG]*m 



(0) 



*(2) 



.(2) 



G * m\;^ + {(%'-Afj,OV^«j' + {5u'-Mu')V,.j}G * m^") + {Sw - Mu')V,> k[dkG] * m 



,(1) 



,(1) 



^(0) 



+iSjy~Mj,,)Vyk[dkG] * m\'> + {5u'~Mu')V,,k{6jj'~Mj,,)V,,i[dkdiG] * m^ 
Introducing the abbreviation Q := (Id — M) V and using the relation (df) * g ~ d{f * g) this can be written as 



(B7a) 
(B7b) 



*(i) 



.(1) 



G*mf> +Qu'd,,{G*m^''^) 



m*f ^ = G * m^f + (Q,j+ Q,,) G * m^") + Qu'd,,{G * mf) + Qjrdj,{G * mf ^) 



'Qjj' di'dj'{G * m 



(0)^ 



Substituting G * m 



(0) 



m' 



(") according to Eq. (B2 1 gives 



(B8a) 
(B8b) 



*(i) 



= G*m\^' +Qu'd,'m*^°'> 



m*f ^ = G * m^) + (Q„- + Q,,) m*(o) + Qwd,, (G * mf'>) + Q ,,,&,. (G * mf^) + i 



Now G * m\ ' 



m, 



<i) 



Qii'di'm*^^^ can be substituted, what leads to 



(B9a) 

■d,,dj-m*^°\ (B9b) 



m*]^' = G * m^^' + (Qij + Qji - Qii'Qjy di-dj')m*'~'^'> + Qti'di'm.y^' + Qjj,dj'm*^^' . 



(BIO) 



Putting these results together and using Eq. M to express m'^'''' in terms of h^*^^ and m'^^ yields the final expressions for the 



noisy moments m 



<k} 



Ko) ^ g, 



,(0) 



*(i) 



G*(/if)m(°)) + Q,,,5,.TO*(") 



f^ = G * ih\fm^'''>) + (Q,j + Q,, - Qu'Qjr d.,dj,)m*<^''^ + Qu'd,,m*/^^ + Q,j,dj,mf'^ 



(Blla) 
(Bllb) 
(Bile) 
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